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,— i Abstract 

A stability analysis is presented of the hydrolysis of methyl isocyanate (MIC) 
using a homogeneous flow reactor paradigm. The results simulate the thermal 
runaway that occurred inside the storage tank of MIC at the Bhopal Union Carbide 
plant in December 1984. The stability properties of the model indicate that the 
thermal runaway may have been due to a large amplitude, hard thermal oscillation 
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initiated at a subcritical Hopf bifurcation. This type of thermal misbehavior cannot 
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be predicted using conventional thermal diagrams, and may be typical of liquid 
thermoreactive systems. 
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^h 1 1. Introduction 

|> 

2 More than 25 years after the Bhopal disaster its horrific legacy is now well- 

3 documented ( |Mishra et al.[ |2009| ), but the causes are still being debated in the 



4 international media. Was the tragedy due to neglect, parsimony, or procrastina- 

5 tion by Union Carbide on safety and maintenance? Ignorance, corruption, sabo- 
e tage and cover-up? Inadequate regulation of urban and industrial development? 
7 Possibly all of the above, but they are putative, secondary or socioeconomic con- 
a tributing factors. (A brief account of the disaster is given in the Appendix.) The 
9 primary cause of the thermal runaway that led to the venting of a poisonous mist 

of methyl isocyanate over the city of Bhopal was physicochemical. In this be- 
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11 lated work I present a simple stability analysis of the thermokinetics of methyl 
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12 isocyanate (MIC) hydrolysis, revealing rogue thermal misbehavior that appears to 

13 be endemic to reactive organic liquids and that cannot be predicted using conven- 

14 tional heat generation/loss rate diagrams. 

15 The Bhopal incident is by no means the only example of disastrous thermal 
runaway occurring in a storage tank. The Seveso accident in Italy in July 1976 
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17 in which large quantities of a toxic dioxin were released into the environment 
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occurred under storage tank conditions (Theophanous 1983 1 , while more or less 



19 minor runaways and explosions due to unforseen reactions in storage tanks and 



vessels are relatively common. The special dangers of chemical storage were 
21 discussed by Gygax ( |1988 ), who pointed out that circumstances favouring heat 



22 accumulation are actually more likely to occur in storage tanks and equipment 

23 parts that are not actively controlled, than in dedicated, well-designed chemical 

24 reactors. 

25 Despite the acknowledged dangers of large-scale storage of reactive chemi- 

26 cals very little has been published on the dynamics of processes that may occur 

27 in storage tanks. Velo et aL] ( 1996[ ) summarize the literature on theoretical and 



28 experimental validations of runaway criteria and parametric sensitivity in batch 

29 reactors and storage tanks. In defining critical conditions they, along with other 

30 authors cited therein, begin with the assumption that storage tanks of small vol- 

31 ume can be modelled as well-stirred batch reactors with linear thermal coupling 

32 to the environment. 

33 However batch reactors have no non-trivial steady states, and there is no gen- 

34 eral theory for determing whether a thermal excursion will grow or decay. Here it 

35 is assumed that the same parameters that govern the stability of a thermoreactive 

36 process in nonequilibrium steady state also govern the stability of thermoreactive 

37 processes in storage tanks. From the comprehensive stability and bifurcation anal- 



ysis of the CSTR that was carried out in Ball ( 1999 ) these parameters are ambient 



temperature, residence time, heat loss, and intensive properties of the reacting sys- 
tem. It is shown in this work that a simple spatially homogeneous, steady state 
approximation can simulate thermoreactive processes in a storage tank with high 
fidelity. 

Another driver for better understanding of thermoreactive processes in liquids 
has emerged recently; this is the use of organic hydroperoxide explosives by ter- 



45 ronsts. 



Table 1: Known relevant data for the MIC thermal runaway at Bhopal. 



Boiling point of MIC at 1 atm 

Density of MIC 

Constant pressure heat capacity of MIC 

Reaction enthalpy for MIC hydrolysis 

Reaction frequency for hydration of MIC 



39.1° C 

0.9599 g/cm 3 at 20° C 

67.7 J/(Kmol) - 1188J/(Kkg) 

-65.1kJ/mol 

3.9 x 10 12 s- : 



(Castro etal. 1985 1 



Activation energy for hydration of MIC 64kJ/mol 



(Castro etal. 1985 1 



Mass of MIC in Tank 610 

Initial temperature inside Tank 610 

Estimated time to criticality 



41 tonnes 
13° C 
4 hours 
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2. MIC chemistry and known relevant data 

A chemical analysis of the residue in the MIC storage tank (Tank 610) sampled 
seventeen days after the event found a variety of condensation products, mainly 
the cyclic trimer ( |D'Silva et al. , 1986). However, experiments to elucidate the 



organic chemistry of the formation of these products indicated that these con- 
densations must have been initiated at temperatures and pressures well above the 
normal boiling point of MIC. Therefore, it is thought that the initial reaction of 
thermokinetic significance was hydration to the unstable N-methyl carbamate, in- 
dicated in grey: 

CH 3 NCO + H 2 ^-l CH3NHCOOH, 

where k(T) is the temperature-dependent (pseudo-)first order rate constant. The 
primary thermal effect that led to the onset of critical conditions is thought to be 
due to the overall reaction 

CH 3 NCO + H 2 ^-l CH3NHCOOH -^ CH 3 NH 2 + C0 2 . 
Some relevant physicochemical data and quantities for MIC hydrolysis in 
Tank 610 are given in table [T] The thermodynamic parameters are taken from 
standard tables. 



3. Thermokinetics of MIC hydrolysis and thermal instability 

The spatially homogeneous flow reactor in which a reactant undergoes a first 
order, exothermic conversion is a simple but elucidatory model for thermoreac- 



tive systems when it is appropriate to ignore convection, because as a dynamical 
system it has non-trivial steady states that can be analysed for stability. The mass 
and enthalpy equations are 

dc 
V— = -Vck(T) + F(c f -c) (1) 

at 

AT 

VC— = (-AH)Vck(T) - (FC + L)(T - T a ). (2) 

at 

V is the reacting volume, c is the reactant concentration, cy is the reactant concen- 
tration in the inflow, F is the volumetric flow rate, C is the averaged volumetric 
specific heat, T is the reaction temperature, AH is the reaction enthalpy, L is the 
linear heat transfer coefficient, T a is the ambient temperature. The temperature- 
dependent (pseudo)-first order reaction rate constant is 

k(T) = A exp(-E/RT), (3) 

where A is the reaction frequency, E is the activation energy, and R is the universal 
gas constant. For numerical and comparative reasons it is useful to work with the 
following dimensionless system corresponding to equations ([l]-[2]), using ([3]): 

d -^ = -xe- lll ' + f{l-x) (4) 

dr 

s d ^ =xe -^_ {ef + £){u _ UaX (5) 

dr 

64 where x = c/c f , r = tA, u = RT/E, f = F/VA, s = CE/c f (-AH)R, £ = 

65 LE/cfVA(-AH)R, u a = RT a /E. Numerical analysis of equations (|4]-|5]) was car- 
ee ried out using rate, fhermochemical and temperature data from tablejTfand values 
67 of the inverse residence time /, heat loss £, and inflow concentration Cf were 
es assigned on the basis of available data. 

69 The reacting mixture self-heats if the rate of reactive heat generation r g ex- 

70 ceeds the linear cooling rate r t . Thermal runaway occurs if r g exceeds r/ beyond 

71 a system- specific threshold; for the hydrolysis of MIC this is taken as the boiling 

72 point of MIC. The steady-state rates from equations (|4]-[5]), are plotted in figure [T] 

73 where the temperature is labeled in dimensional units. The reaction self- heats 

74 until the reactor temperature T reaches the steady state temperature of ~305 K at 

75 which the heating and cooling rates are balanced. Since the boiling point of MIC 

76 is 312 K, on the basis of this diagram we would not expect a thermal runaway 

77 to develop, even when the ambient temperature is allowed to drift slowly up to 

78 292 K. 
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Figure 1: Rates of reactive heat generation r g (red) and heat loss r; (blue) versus T from equations 
(ffl-B). u a = 0.0379 (corresponding to T a = 292 K), f =1.1,1 = 700, s = 10. 
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However thermal balance diagrams such as that in figure[T]can be dangerously 
misleading because they infer stability rather than assess stability rigorously, al- 
though such diagrams are often used in chemical reactor engineering. The steady 
states, periodic solutions, and stability analysis of equations (j4]-[5j) were computed 
numerically (Doedel 2010[ ) and yielded a dramatically different picture of the the 
thermal stability of MIC hydrolysis. Figure [2] is a bifurcation diagram in which 
the steady states and the amplitude envelope of periodic solutions are plotted as 
a function of T a . The steady state is stable at T a « 2S6K, the temperature at 
which the tank of MIC had been held for several months. As T a is quasistatically 
increased the reaction temperature T increases slowly, but at T a = 290.15^ the 
stability analysis flags an abrupt change in the nature of the solutions. At this point 
the steady state solutions lose stability to a Hopf bifurcation and the hard onset 
of a high amplitude thermal oscillation ensues. Clearly, at T a = 292 K we have 
catastrophic thermal runaway, contrary to the inadequate prediction given by fig- 
ure^ (In the resulting superheated liquid the exothermic condensation reactions 
would increase the temperature even further.) 

This is quite different from classical ignition of a thermoreactive system, which 
occurs at a steady-state turning point. The dynamics of oscillatory thermal run- 
away can be understood by studying a close-up of the region around the Hopf 
bifurcation Hi in figure [2j This is shown in figure [3] Hi is subcritical and the 
emergent limit cycle is unstable. The amplitude envelope of the unstable limit cy- 
cles is marked with a thin dotted line; they grow as T a is decreased. At the turning 
point of the periodic solution branch marked with a large asterisk the limit cycles 
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Figure 2: Bifurcation diagram. Stable steady states are plotted with solid blue line, unstable steady 
states with dashed red line, and the amplitude envelope of periodic solutions is marked with thin 
dotted magenta line. H\ and H^ label the Hopf bifurcation and the large # marks the change in 
stability of the limit cycles. / = 1.7, ( - 700, s - 10. 



become stable. Thermal runaway may occur if there are significant perturbations 
while T a is within the regime of unstable limit cycles, and it must occur when 
T a drifts above H\ . The arrow indicates the rapid thermal excursion, in principle 
to the stable limit cycle but in reality the reactant has vaporised, the pressure has 
soared beyond the safety limits of the tank, and the system must vent since the 
peak temperature is far above the boiling point of MIC. 
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4. Discussion 

Although in principle classical thermal ignition is possible for MIC hydrolysis, 
in practice the presence of oscillatory instability is all-pervasive and dominant in 
this system. This can be appreciated by inspection of figure |4} a plot of the loci of 
the steady state turning points and the Hopf bifurcations of equations (|4]-[5]) over 
the two parameters u a and the inverse residence time /. In the filled region thermal 
runaway will always be oscillatory. The bistable regime occurs at very high flow 
rates (short residence times). Here classical ignition at a steady state turning point 
to a stable steady state may occur, but the oscillatory instability is still present and 
oscillatory thermal runaway occurs (perversely) as the ambient temperature u a is 
reduced. 

The tendency for oscillatory, non-classical thermal runaway may be typical of 
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Figure 3: Close-up of the region around the Hopf bifurcation Hi in figure [2] 
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low-boiling, exothermically reactive organic liquids. In the work of Ball and Gray 



(1995) the hydration of 2,3-epoxy-l-propanol in a CSTR was found to exhibit 
similar non-classical thermal misbehavior. The design of safe storage systems 
for such liquids should focus on damping the oscillatory instability, rather than 
shifting a classical ignition point using over- simplistic heat generation/loss rate 
diagrams like figure [TJ 

In recent years low-boiling, exothermically reactive liquids such as a variety 
of organic peroxides have been used as explosives by terrorists, and their potential 
use as murder weapons on aircraft is the reason for current restrictions on the 
liquids that passengers may carry on-board. It is possible that peroxide explosions 
are due to oscillatory instability rather than classical thermal ignition. Figure [5] 
shows the two-parameter bifurcation diagram using thermokinetic parameters for 
the decomposition of cumene hydroperoxide (Wu and Shu, 2010[ ) in equations 
(|4]-[5]). In understanding this type of explosion and for deactivating improvised 
explosive devices that employ such liquids it is clearly necessary to understand 
oscillatory thermal instability. 

Is it realistic to model a reacting volume inside Tank 610 as a well-stirred flow 
reactor? Yes, on a timescale over which the reacting volume remains relatively 
constant and gradientless relative to the much faster rate of reaction. Consider a 
volume V of liquid in the tank, near the source of the water ingress, within which 
the hydration reaction takes place. The reacting volume V grows as water flows 
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Figure 4: The locus of Hopf bifurcations is marked with a solid line, the locus of steady-state 
turning points is marked with dashed line. I - 700, s — 10. 
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into it and reaction proceeds. This volume can be regarded as the "reactor", while 
the uncontaminated, non-reacting MIC in the rest of the tank can be regarded as 
the "coolant" to which heat is transferred linearly. The reacting volume expands 
until the rate of inflow matches the rate of outflow of reactant and products from 
the reacting volume. 

In other words, for the purposes of this model in which the focus is on the 
dynamics we circumscribe a volume in which the spatial gradients are insignifi- 
cant in comparison to the time evolution, and therefore can be neglected. If this 
approximation does not hold, then we are free to reduce the circumscribed vol- 
ume until it does. There is nothing particularly artificial or manipulative in doing 
this; it is just a simplest case scenario for which the powerful tools of stability and 
bifurcation theory can be applied. 

The Union Carbide plant at Bhopal has been derelict since the disaster and 
MIC is no longer in use as a bulk chemical, but we cannot afford to close the 
book on potential problems in storage tanks. In chemical plants and storage sites 
around the world there are tanks containing reactive organic liquids that have sim- 
ilar physical properties to MIC and undergo reactions with similar thermokinetics. 
Many other chemicals are stored as bulk commodities that can undergo thermally 
unstable polymerization, oxidation or decomposition reactions. Analysis of sim- 
ple gradientless models of processes that may occur in storage tanks provides 
insight into the physical basis of rogue thermal behavior and a starting point for 
improved design of safe storage systems. Concomitantly there is a need for more, 
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Figure 5: The 2-parameter bifurcation diagram from equations fflpl) applied to the decompo- 
sition of cumene hydroperoxide also has a large region within which thermal runaway is due to 
oscillatory instability. The locus of Hopf bifurcations is marked with a solid line, the locus of 
steady-state turning points is marked with dashed line. I - 700, s - 20. 



163 and more accurate, thermokinetic data to use in such models. It is notable that the 

164 data for MIC hydration/hydrolysis in table[T]were obtained as part of an unrelated 

165 research project initiated before the Bhopal disaster. 



166 5. Conclusions 

167 1. A belated investigation of the thermal runaway in Tank 610 that led to the 
lea Bhopal disaster has been carried out by modelling the hydrolysis of MIC in 

169 the well-stirred limit and analysing the stability of solutions of the dynami- 

170 cal system model. 

171 2. Thermal runaway occurs due to the onset of a hard thermal oscillation at 

172 a subcritical Hopf bifurcation. Classical thermal ignition at a steady state 

173 turning point may occur in principle but over the thermal regime of interest 

174 the system is dominated by oscillatory instability. 

175 3. This non-classical oscillatory thermal misbehavior may be generic in reac- 

176 tive organic liquids. The results yield valuable intelligence about the causes 

177 of themal runaway that may inform improved designs of storage systems 

178 for large quantities of toxic and reactive substances. 

179 4. These results may also inform better management of organic peroxide based 

180 explosives. 



181 Appendix 

182 The following brief account of the Bhopal disaster has been compiled from the 
following sources: |Forman| ( |1985[ ), [Lepkowski| ( |1994[ ), [Shrivastava| ( |1987[ ), [Weir 
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184 ( fT987| ) and [Abbasi and Abbasil ( |2005] ). 



185 The Union Carbide plant at Bhopal carried out the production of carbaryl, 

186 an agricultural insecticide that has been used widely throughout the world since 

187 1945. Methyl isocyanate, a low-boiling, highly reactive and extremely toxic liquid 
lea used in the synthesis of carbaryl, was stored in an underground stainless steel tank 

189 (Tank 610) which was encased in a concrete shell. The temperature of the 41 

190 tonnes of MIC in Tank 610 was 12-14° C rather than the recommended 0-4° C 

191 because the refrigeration unit had been non-operational for several months. On 

192 the evening of December 2 1984 a worker had been sent to hose out a nearby 

193 tank. The hose was left running unattended, and it is believed that a faulty valve 

194 allowed entry of water into the connected Tank 610. (Union Carbide disputes this, 

195 asserts that nothing was wrong with its equipment and procedures, and argues 

196 that sabotage by a disaffected employee must have caused the disaster.) By 1 1:30 

197 pm, when workers detected lachrymose whiffs of leaking MIC, water had been 

198 running into Tank 610 for at least four hours. Although a slow rise in temperature 

199 and pressure in the tank had been noted, the early signs of trouble were not acted 

200 upon. Shortly after 1 1:30 pm the contents of the tank reached thermal criticality 
and began escaping as vapor from the flare tower. 

Downwind of the flare tower lay the crowded suburbs and shantytowns. Most 
of the fluid in the tank streamed from the tower then drifted low over the city and 
sank and seeped in deathly mist in lungs in eyes, a period to sleep and swift ar- 

205 rest of retreat. More than 3000 lives were claimed immediately and many tens of 

206 thousands through the subsequent days and months and years lost their lives or 
their health to the poison's effects, and the dead are still being counted. 
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